;; Tests of the non-robust computational geometry functions
;; Copyright (C) 2022  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.delaunay.bowyer-watson-test
  (:require
   [clojure.test :refer [deftest is testing]]
   [mazama.compgeom.convex-hull.monotone-chain :refer [convex-hull]]
   [mazama.compgeom.core :refer [incircle point-orient2d]]
   [mazama.compgeom.delaunay.bowyer-watson :refer [delaunay ghostbuster triangles]]))

(defn- read-number [s]
  (clojure.edn/read-string s))

(defn- read-node-line [line]
  (let [[_ x y] (clojure.string/split line #"\s+")]
    (vector (read-number x) (read-number y))))

(defn- ignore? [line]
  (or (empty? line)
      (clojure.string/starts-with? line "#")))

(defn- normalize [file]
  (if (.endsWith file ".gz")
    (-> file
        clojure.java.io/input-stream
        java.util.zip.GZIPInputStream.)
    file))

(defn read-node-file
  "Read a node file.
  https://www.cs.cmu.edu/~quake/triangle.node.html"
  [file]
  (with-open [rdr (-> file normalize clojure.java.io/reader)]
    (->> (line-seq rdr)
         (map clojure.string/triml)
         (filter (complement ignore?))
         (drop 1)
         (map read-node-line)
         vec)))

(defn- read-ele-line [line]
  (let [[_ a b c] (clojure.string/split line #"\s+")]
    (vector (read-number a) (read-number b) (read-number c))))

(defn read-ele-file
  "Read an ele file.
  https://www.cs.cmu.edu/~quake/triangle.ele.html"
  [file]
  (with-open [rdr (clojure.java.io/reader file)]
    (->> (line-seq rdr)
         (map clojure.string/triml)
         (filter (complement ignore?))
         (drop 1)
         (map read-ele-line)
         vec)))

(defn- rotate
  "Rotate the triple so that the smallest value comes first."
  [[a b c]]
  (let [m (min a b c)
        i (cond
            (== m a) 0
            (== m b) 1
            :else 2)
        t [a b c]]
    [(t i)
     (t (mod (+ i 1) 3))
     (t (mod (+ i 2) 3))]))

(defn convert-mesh-to-ele-format [mesh S]
  (let [P (into {} (map vector S (iterate inc 1)))]
    (map (fn [[a b c]] [(P a) (P b) (P c)]) mesh)))

(defn print-mesh-using-ele-format [mesh S]
  (let [Q (convert-mesh-to-ele-format mesh S)
        I (->> Q
               (map rotate)             ;optional
               sort                     ;optional
               (map cons (iterate inc 1)))]
    (println (count I) 3 0)
    (doseq [[i a b c] I] (println i a b c))))

(defn write-ele-file
  "Write an ele file."
  [file mesh S]
  (with-open [wtr (clojure.java.io/writer file)]
    (binding [*out* wtr]
      (print-mesh-using-ele-format mesh S))))

(defn- rationalify [S]
  (map #(mapv rationalize %) S))

(defn- euler
  "Euler's formula for the number of triangles in a Delaunay triangulation.
  Page 76 https://www.cs.umd.edu/class/spring2020/cmsc754/Lects/cmsc754-spring2020-lects.pdf"
  [S]
  (- (* 2 (count (set S)))
     2
     (count (convex-hull S true))))    ;true -> all boundary points are on the hull

(defn- circumcircle-intersections [T S]
  (for [t T s S :when (pos? (incircle t s))] [t s]))

(defn- check-delaunay
  "Is `mesh` a Delaunay triangulation of the points in `S`?
  `S` is assumed to have at least three non-colinear points."
  [mesh S]
  (let [T (-> mesh ghostbuster triangles)
        S (set S)]
    ;; T and S have the same points
    (is (= (sort S)
           (->> T (reduce into #{}) sort)))
    ;; Check the number of triangles
    (is (== (count T) (euler S)))
    ;; Check every triangle is counterclockwise
    (is (every? #(let [[p q r] %] (pos? (point-orient2d p q r))) T))
    ;; Check for no circumcircle intersections with points in S
    (is (empty? (circumcircle-intersections T S)))))

(deftest not-delaunay
  ; Testing check-delaunay
  (testing "circumcircle intersections"
    (let [S [[0 0] [4 0] [3 -3] [2 3]]
          T [[[2 3] [3 -3] [0 0]] [[[4 0] [3 -3] [2 3]]]]]
      (is (not (empty? (circumcircle-intersections T S)))))))

(deftest degenerate-test
  (testing "nil"
    (is (empty? (delaunay nil))))
  (testing "empty"
    (is (empty? (delaunay []))))
  (testing "one vertex"
    (let [S [[332 0]]]
      (is (empty? (delaunay S)))))
  (testing "two vertices"
    (let [S [[3 -4] [-37 8]]]
      (is (empty? (delaunay S)))))
  (testing "three colinear vertices"
    (let [S [[1 2] [3 4] [5 6]]]
      (is (empty? (delaunay S))))))

(deftest duplicate-vertices-test
  (testing "duplicate vertex among the first triangle"
    (let [S [[0 0] [1 1] [1 0] [-2 -1] [1 1]]]
      (check-delaunay (delaunay S) S)))
  (testing "duplicate vertex not among the first triangle"
    (let [S [[0 0] [1 1] [1 0] [2 3] [7 3] [-1 2] [2 3]]]
      (check-delaunay (delaunay S) S))))

(deftest triangle-test
  ; From http://people.eecs.berkeley.edu/~jrs/274f09/proj.html
  (let [tests ["http://people.eecs.berkeley.edu/~jrs/input/tri.node"
               "http://people.eecs.berkeley.edu/~jrs/input/4.node"
               "http://people.eecs.berkeley.edu/~jrs/input/box.node"
               "http://people.eecs.berkeley.edu/~jrs/input/spiral.node"
               "http://people.eecs.berkeley.edu/~jrs/input/flag.node"
               ; The following fails. Not robust enough to handle floating point. Rational version below...
               #_"http://people.eecs.berkeley.edu/~jrs/input/grid.node"
               "http://people.eecs.berkeley.edu/~jrs/input/dots.node"
               "http://people.eecs.berkeley.edu/~jrs/input/ladder.node"
               "http://people.eecs.berkeley.edu/~jrs/input/633.node"
               ; The following is very slow (because check-delaunay is very slow).
               ; It takes about 6 minutes but it's worth running every once in a while.
               #_"http://people.eecs.berkeley.edu/~jrs/input/ttimeu10000.node.gz"]]
    (doseq [t tests] (let [S (read-node-file t)]
                       (check-delaunay (delaunay S) S)))))

(deftest rational-triangle-test
  ; Triangle tests that need to be rationalized
  (let [tests ["http://people.eecs.berkeley.edu/~jrs/input/grid.node"]]
    (doseq [t tests] (let [S (rationalify (read-node-file t))]
                       (check-delaunay (delaunay S) S)))))

(deftest bug-test
  (testing "Bug introduced by commit c0e78ea65b9739322fd200a495609b483852dbc6"
    ; The :mru was left in an incorrect location (at the end of the initial colinear points).
    (let [S [[0 0] [0 1] [0 2] [1 2] [2 0]]]
      (check-delaunay (delaunay S) S))))
